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PREFACE 


This  Note  covers  work  done  on  the  problem  of  modeling  and 
estimating  the  demand  for  aircraft  spare  parts  while  the  author  was  a 
consultant  to  The  Rand  Corporation  in  the  summer  of  1983.  The  research 
began  under  the  Project  AIR  FORCE  Resource  Management  Program  study 
entitled  "The  Driving  Inputs  and  Assumptions  of  Stockage  Assessment 
Models."  The  Note  is  being  published  as  part  of  a  follow-on  study 
entitled  "Enhancing  the  Integration  and  Responsiveness  of  the  Support 
System  to  Meet  Wartime  and  Peacetime  Uncertainties." 

This  work  grew  out  of  a  study  of  parts  failure  data  from  several 
Air  Force  units,  strongly  indicating  that  either  current  models  of  part 
failure  behavior  or  prevailing  beliefs  about  the  inherent  stability  of 
this  behavior,  or  both,  were  wrong.  The  results  here  support  this 
indication,  but  they  are  of  general  interest  for  modeling  and  estimating 
nonhomogeneous  Poisson  processes. 

This  Note  will  be  of  interest  to  those  concerned  with  forecasting 
inventory  requirements  in  the  Department  of  Defense  or  industry. 
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SUMMARY 


Mathematical  models  are  commonly  used  to  study  tie  performance  of 
the  Air  Force's  spare  parts  supply  and  repair  systems.  But  accurate 
evaluations  of  supply  policies  are  not  possible  without  accurate  models 
of  the  supply  system,  and  models  that  understate  the  variability  in  the 
supply  system  will  bias  evaluations  in  favor  of  policies  that  rely  on 
accurate  predictions  of  part  failures.  This  Note  examines  the  model  for 
part  failures  used  in  The  Rand  Corporation  Supply  System  model 
Dyna-METRIC. 

Section  I  gives  a  short  description  of  Dyna-METRIC  and  then  a  longer 
description  of  common  probability  models  for  part  failures.  The 
strengths  and  weaknesses  of  these  models  are  discussed,  with  particular 
attention  to  sources  of  variability  in  observed  part  failure  behavior 
that  the  models  do  not  appear  to  capture. 

Section  II  has  two  purposes.  The  first  is  to  examine  the  plausibility 
of  Dyna-METRIC* s  current  probability  model  for  part  failures  in  the 
light  of  some  new  Air  Force  data.  This  model  treats  the  number  of 
failures  of  a  particular  part,  at  a  particular  air  base,  in  a  time 
period  of  given  length,  as  a  Poisson  random  variable  with  mean  Xf,  where 
f  is  the  total  number  of  hours  flown  in  that  time  period,  and  X  is  an 
unknown  constant  characteristic  of  the  part  and  airbase.  The  second 
purpose  is  to  derive  some  new  properties  of  V,  a  commonly  used  estimator 
of  the  variance-to-mean  ratio,  under  thre_  different  probability  models 
for  part  failures.  These  properties  of  V  indicate  strongly  that 
Dyna-METRIC1 s  current  probability  model  does  not  permit  enough 
variability  to  credibly  explain  the  Air  Force  data.  Further,  although 
the  data  do  indicate  that  it  is  preferable  to  model  mean  part  failures 
as  a  +  Bf  for  c  and  B  unknown  constants  and  a  not  necessarily  zero,  a 
Poisson  model  with  this  mean  still  does  not  allow  enough  variability.  A 
negative  binomial  model  with  mean  a  +  Bf  is  preferred. 

Finally,  Sec.  II  shows  that  V  is  always  biased  low  for  the 
probability  model  in  which  it  is  intended  to  be  used — i.e.,  with  part 
failures  distributed  as  negative  binomial  random  variables  with  mean  Xf 

PREVIOUS  PACE 
»s  blank 


and  variance-to-mean  ratio  p.  This  bias  is  an  increasing  function  of  p 
for  a  fixed  number  of  total  expected  failures  and  can  be  very  large  for 
large  p. 

Section  III  contains  suggestions  for  estimating  the  parameters  of  the 
models  recommended  in  Sec.  II.  Maximum  likelihood  estimates  are  suggested 
because  they  are  tractable  and  because  they  appear  to  solve  the  bias 
problem  noted  above. 
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I.  THE  REAL  WORLD  PROBLEM,  SOME  CURRENT  MODELS 
AND  SOME  MODELING  ISSUES 


INTRODUCTION 

One  of  the  missions  assigned  to  the  Air  Force  is  to  be  able  to 
relocate  fighting  units  and  associated  service  and  supply  facilities  to 
far-flung  places  on  short  notice,  to  meet  a  variety  of  policy 
commitments  -  A  remarkable  amount  of  materiel  is  needed  to  keep  such 
units  operational,  and  airlifting  capacity  is  limited  and  expensive. 
Because  of  this  limited  capacity,  several  policy  choices  must  be  made. 
One  such  choice  is  how  best  to  handle  spare  parts  supply  requirements — 
i.e.,  is  it  possible  to  design  spare  parts  kits  that  don't  take  up  toe 
much  airlift  space  but  are  sufficient  to  maintain  an  adequate  number  of 
operational  planes  for  some  desired  length  of  time;  or  would  it  be 
preferable  to  have  a  "responsive"  supply  system  in  which,  for  example,  a 
dedicated  fleet  of  planes  hauls  spare  parts  around  as  needed,  reducing 
spare  parts  kits  airlifted  with  the  units  accordingly?  Given  the  choice 
of  a  general  approach  to  supplying  parts,  how  often  should  shipments  be 
made,  how  large  should  these  shipments  be?  And  so  on. 

These  kinds  of  questions  have  long  been  examined  with  the  help  of 
mathematical  models  intended  to  simulate  the  wartime  performance  of  the 
spare  parts  supply  and  repair  system.  One  model  developed  at  Rand  and 
in  wide  use  within  the  Air  Force  is  called  Dyaa-METRIC  (for  Dynamic 
Multi-Echelon  Technique  for  Recoverable  Item  Control);  henceforth  I  will 
often  refer  specifically  to  Dyna-METRIC,  but  its  features  are  not 
atypical  of  such  models. 

These  models  attempt  to  characterize  the  (1)  numbers  of  failures  of 
parts  on  planes ,  where  proneness  to  failure  is  allowed  to  depend  on  the 
intensity  of  use  as  dictated  by  a  war  scenario  (thus  "Dynamic”);  (2) 
repair  times  of  those  parts  that  can  be  repaired,  either  at  the  forward 
airbase  or  at  some  more  centralized  rearward  repair  facility  (thus 
"multi-echelon");  (3)  transportation  times  and  other  delays ,  where 
relevant;  and  (4)  stocking  requirements ,  the  levels  of  stocks  of  parts 
kept  at  the  forward  base  and  at  the  more  centralized  repair  facility, 
and  levels  of  and  timing  of  supplies  of  new  parts. 


The  general  scheme  of  the  Dyna-METRIC  model  is  depicted  in  Fig.  1, 
which  is  a  simplified  version  of  Fig.  ]  in  Hillestad  (1982).  In 
Dyna-METRIC,  the  numbers  of  part  failures  are  modeled  as  random 
variables  whose  probabi lity  distributions  depend  on  the  number  of  hours 
flown,  the  measure  of  intensity  of  use.  To  be  more  specific,  consider  a 
particular  Air  Force  unit,  and  a  part  i,  which  is  in  all  of  the  planes 
in  that  unit.  Let  X^(i)  denote  the  number  of  failures  on  part  i  in  that 
un^  *n  time  period  j.  In  Dyna-METRIC,  the  parameters  of  the 
probability  distribution  of  X^(i)  are  assumed  to  depend  on  total  hours 
flown  by  the  unit  in  the  jth  time  period.  The  effects  of  different  w’ar 
scenarios  on  failures  are  examined  by  manipulating  the  number  of  hours 
flown,  the  number  of  planes  used,  and  other  factors  affecting  the 
distributions  of  X^(i)  for  the  various  parts. 


Fig.  1  —  Simplified  scheme  of  parts  repair  and  supply  model 
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Repair  times  for  individual  failed  parts  are  treated  as  random 
variables  having  exponential  or  degenerate  distributions,  and  assumptions 
are  made  about  the  Drobabilities  that  parts  can  be  repaired  only  at  the 
rearward  centralize!  facility.  The  exponential  distribution  was  chosen 
for  mathematical  convenience,  but  it  seems  to  be  a  satisfactory 
assumption.  Both  under  the  assumption  of  an  infinite  number  of  repair 
stands  (Crawford,  1981)  and,  as  far  as  is  known,  under  the  assumption  of 
a  finite  number  of  repair  stands,  the  supply  system  performance 
measurements  produced  by  Dyna-METRIC  are  fairly  insensitive  to  the 
distributional  forms  assumed  for  repair  times. 

The  effects  of  various  stockage  policies  can  be  introduced  by 
vs  ying  purchases  from  vendors  and  decisions  about  allocating  spare 
parts  to  the  stock  of  the  forward  bases  or  the  rearward  facility.  It  is 
assumed  that  stocking  policy  is  nonstochastic,  and  these  policies  are 
represented  by  changes  at  prescribed  times  in  the  numbers  in  stock  at 
the  bases  and  the  rearward  facility. 

The  purpose  of  these  models  is  to  connect  dollars  and  cents  policy 
decisions  to  aircraft  performance  of  assigned  missions.  Dyna-METRIC 
and  other  models  produce  several  performance  measures.  Some  are 
inventory  measures,  such  as  the  probability  distribution  of  the 
backorder  across  all  parts,  the  probability  that  a  failure  of  a 
particular  part  will  be  filled  immediately  from  stock,  and  so  on.  These 
measures  are  derived  analytically,  based  on  the  distributions  assumed 
for  failures  and  repair  times.  More  relevant  to  the  policy  questions, 
if  mission  requirements  can  be  specified  in  terms  of  numbers  of  planes 
with  particular  capabilities,  it  is  possible  to  derive  a  probability 
distribution  for  mission  capability  from  the  already  noted  inventory 
measures . 

MODELING  FAILURES  (DEMANDS) 

These  performance  measures  are  no  better  than  the  model  that 
generates  them.  The  particular  area  of  concern  here  is  modeling  the 
numbers  of  failures,  or  demands,  and  estimating  the  parameters  in  those 
models. 


A  common  assumption  is  that  the  number  of  failures  X  in  a  period  of 
given  length  has  a  Poisson  distribution  with  parameter  y,  in  which  case 


P(X  =  k)  =  e'VVk!  for  k  =  0,  1,  2,  - 

Under  this  assumption,  both  the  mean  and  variance  of  X  are  equal  to  y, 
so  that  the  variance-to-mean  ratio  p  =  var(X)/E(X)  is  equal  to  1-  To 
model  applications  in  which  p  exceeds  one,  a  negative  binomial 
distribution  is  often  hypothesized.  Here, 

P(X  =  k)  =  (k  +kr  X) ( 1  -  p)kpr  for  k  =  0,  1,  2,  ..., 

where  r  >  0,  0  <  p  <  1.  In  this  case,  E(X)  =  r(l  -  p)/p  and  var(X)  = 

2 

r(l  -  p)/p  ,  so  that  p  =  1/p  >1.  As  is  well  known,  the  Poisson 
distribution  is  a  limit  of  negative  binomial  distributions  as  p  ■*  1  and 
r  ■*  «•  in  such  a  way  that  r(l  -  p)  remains  fixed  at  y  (Feller,  1968,  p. 
172). 

The  current  version  of  Dyna -METRIC  assumes  that  demand  for  a 
particular  part  for  a  particular  Air  Force  unit  in  some  period  of  given 
length  is  either  Poisson  or  negative  binomial.  When  the  Poisson 
distribution  Is  used,  its  mean  is  assumed  to  equal  Xf,  where  X  is  an 

unknown  constant  peculiar  to  the  part  and  unit,  and  f  is  the  number  of 

hours  flown  by  the  unit  in  the  given  period.  When  the  negative  binomial 
distribution  is  used,  its  mean  is  also  assumed  to  be  Xf,  and  its 

variance-to-mean  ratio  is  assumed  to  be  p,  so  that  r  =  Xf/(p  -  1)  and 

P  =  1/P- 

The  negative  binomial  is  a  "compound  Poisson"  distribution  in  both 
senses  in  which  that  term  is  cossnonly  used.  In  the  first  sense,  usually 
associated  with  probabilists ,  the  random  variable  =  X^  +  X,  +  —  + 
Xj.  can  be  shown  to  be  negative  binomial  if  the  X^  are  independently  and 
identically  distributed  with 

P(Xj  =  n)  =  (1  -  p)n/(-n  £n  p),  n  =  1,  2,  ..., 

for  0  <  p  <  1,  and  N  is  independently  distributed  as  a  Poisson  random 
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variable  with  y  =  -  r  in  p  (Feller,  1968,  pp.  268-269  and  286-291).  In 
the  second  sense,  usually  associated  with  statisticians,  if  the 
conditional  distribution  of  a  random  variable  X,  given  its  mean  y,  is 
Poisson,  and  y  itself  has  a  gamma  distribution  with  density 

g(y)  =  &"V"V'j/B/r(a), 

then  unconditionally  X  has  a  negative  binomial  distribution  with  r  =  o 
and  p  =  (8  +  1)  1  (Hogg  and  Craig,  1970,  pp.  99-211).  More  generally, 
compound  Poisson  distributions  are  defined  by  replacing  P(X^  =  n)  with 
an  appropriate  distribution  in  the  first  sense,  or  by  replacing  the 
gamma  distribution  appropriately  in  the  second  sense. 

Aside  from  mathematical  tractability,  the  use  of  compound  Poisson 

models  has  been  rationalized  by  at  least  two  substantive  considerations. 

First,  it  follows  from  Levy's  work  on  infinitely  divisible  distributions 

that  any  arrival  process  where  the  numbers  of  arrivals  in  disjoint  time 

intervals  are  independent  is  Poisson  or  compound  Poisson  in  the  first 

sense  above  (Feller,  1968,  pp.  289-290;  Crawford,  1981,  p.  10).  Thus  an 

assumption  of  independence  of  failures  in  disjoint  time  intervals 

justifies  using  compound  Poisson  models,  although  it  does  not  justify 

the  use  of  any  particular  compound  Poisson  model.  Second,  under  some 

conditions  the  Poisson  distribution  is  a  good  approximation  to  the 

distribution  of  the  s,,-»  P  =  X^  +  —  +  ^  of  n  mutually  independent 

with  distributions  P(X^  =  1)  = 

If  the  probabilities  p^  depend  on  n  in  such  a  way  that  the  largest  p^ 

tends  to  zero  but  the  sum  p,  +■  —  +  p  =  X  remains  constant,  then  in 

*n 

the  limit  as  n  +  «,  has  a  Poisson  distribution  with  parameter  X 
(Feller,  1968,  p.  282).  Thus  for  large  n  and  moderate  values  of  X,  the 
distribution  of  can  be  approximated  by  a  Poisson  distribution.  This 
justification  of  a  Poisson  model  has  appeal  for  modeling  failures  of 
such  parts  as  landing  gear,  for  which  it  is  more  natural  to  assume  that 
a  failure  will  occur  with  seme  fixed  probability  on  each  mission  than  to 
assume  a  Poisson  failure  distribution  for  landing  gear  directly. 

Both  of  these  substantive  considerations  rely  on  an  assumption  of 
independence.  In  some  situations  this  may  be  clearly  inappropriate,  for 
example,  for  parts  that  age  rapidly  or  have  very  high  failure  rates 


random  variables 


=  1  -  P(Xk  “  0)- 


-  6  - 


(Crawford,  1981,  p.  11).  However,  loosening  the  independence  assumption 
would  require  extensive  reworking  of  the  probability  theory  underlying 
the  failure  models,  and  currently  the  indeper ience  assumption  is  not 
believed  to  be  inaccurate  enough  to  justify  t.iis  reworking.  Henceforth, 
we  will  assume  that  failures  follow  a  compound  Poisson  distribution. 

At  this  point,  several  modeling  issues  can  be  raised: 

1.  As  already  noted,  independence  of  numbers  of  failures  in 
disjoint  time  intervals  is  assumed.  This  will  not  be  pursued 
further . 

2.  Is  the  failure  distribution  really  stationary  in  peacetime? 

For  example,  is  there  seasonal  variation  in  the  failure 
behavior  of  certain  kinds  of  parts?  Available  data  indicate 
that  there  is  some  seasonal  variation,  and  it  is  not  difficult 
to  imagine  that  there  would  be  some  ix  extreme  climates - 
Stationarity  will  henceforth  be  assumed. 

3.  Even  assuming  independence,  is  there  any  reason  other  than 
convenience  to  assume  that  the  appropriate  compound  Poisson  is 
a  negative  binomial  and  not  some  othei  compound  Poisson?  There 
has  been  some  investigation  indicating  that,  for  real  aircraft 
failure  data,  the  variance-to-mean  ratio  increases  with  the 
mean,  which  is  not  a  property  of  the  negative  binomial  although 
it  is  a  property  of  other  compound  Poisson  distributions.  The 
actual  import  of  this  apparent  finding  is  unclear,  because  it 
may  simply  be  an  artifact  of  the  estimator  used. 

4.  The  relationship  between  mean  demands  and  flying  hours  is  not 
well  understood.  First,  it  is  not  at  all  clear  that  flying 
hours  are  an  appropriate  "clock" — e.g. ,  for  landing  gear  what 
matters  is  not  how  long  the  plane  is  i*i  the  air,  but  how  often 
it  lands;  also,  some  radar  parts  spend  substantial  amounts  of 
time  switched  on  and  running  while  the  plane  is  on  the  ground, 
so  that  flying  hours  understate  the  actual  intensity  of  use- 
Second,  there  is  no  particular  reason  to  assume  that  the 
relationship  between  flying  hours  and  mean  demands  goes  through 
the  origin.  If  the  planes  simply  sat  in  the  hangars,  some 
failures  would  occur  anyway.  Also,  if  this  relationship  were 
nonlinear  with  a  positive  second  derivative,  and  it  was  desired 
to  approximate  the  relationship  for  higher  numbers  of  flying 
hours,  a  linear  approximation  might  be  appropriate,  but  its 
intercept  would  be  negative-  Third,  there  is  no  particular 
reason  to  assume  that  mean  demands  are  a  linear  nondecreasing 
function  of  flying  hours.  On  the  contrary,  there  is  evidence 
that  high  sortie  rates  can  actually  improve  the  "health"  of 
some  aircraft  parts.  A  study  by  Crawford  and  Kamins 
(forthcoming)  of  objective  measures  of  the  health  of  components 
in  the  fire  control  and  weapons  delivery  system  of  F-I6s, 
gathered  during  a  surge  of  aircraft  activity  during  an 
exercise,  found  that  reported  rates  of  malfunctions,  as 


measured  by  the  aircraft  systems,  actually  declined  comoared 
with  rates  reported  in  periods  of  lower  activity  before  and 
after  the  surge. 


5.  The  available  data  are  somewhat  less  than  ideal.  The  Air  Force 
Logistics  Command  accumulates  failure  data  for  each  part  over 
all  bases  worldwide,  and  disaggregated  data  are  not  currently 
available,  as  far  as  I  know.  The  base  data  available  for  this 
study  were  collected  at  Rand  for  20  parts  with  fairly  high 
demand  rates,  on  five  airbases,  over  at  most  13  quarters. 

Among  the  problems  with  these  data  were:  (a)  "managed  demand," 
cr  possible  changes  in  reported  numbers  of  failures  associated 
with  management  decisions  related  to  reporting  or  servicing; 
and  (b)  a  small  number  of  observations  for  each  part/base 
combination. 

6.  There  was  a  manifest  lack  of  enthusiasm  among  data  collectors 
(mechanics),  the  utility  of  the  data  not  having  been  strongly 
impressed  on  them. 


This  Note  will  consider  only  the  third  of  these  issues,  the  nature 
of  the  compounding  distribution,  and  the  fourth,  the  relationship 
between  flying  hours  and  the  mean  numbers  of  failures.  This  emphasis 
does  not  imply  a  judgment  that  the  other  issues  are  unimportant.  In 
fact,  they  are  so  poorly  understood  that  their  importance  cannot  be 
assessed  satisfactorily.  But  in  the  consideration  of  competing  spare 
parts  supply  systems,  it  is  essential  to  understand  how  well  we  can 
predict  levels  of  parts  failures  under  wartime  conditions;  issues  that 
are  not  well  understood  cannot  be  dismissed. 


II-  SOME  NOTES  ON  AN  ESTIMATOR  OF 
VARIANCE-TO-MEAN  RATIO 


INTRODUCTION 

This  section  describes  some  properties  of  an  estimator  of  variance- 
to-mean  ratio  called  V.  Let 

=  number  of  demands  in  period  i  (i=  1,  2,  . ..,  n), 

N  =  LXi  =  total  number  of  demands  over  all  n  periods, 
f^  =  number  of  flying  hours  in  period  i  (i  =  1,  2,  n), 

T  =  Ef^  =  total  number  of  flying  hours  in  all  periods. 

If  N  >  0,  V  is  defined  to  be 

V  =  s2/x, 

where  X  =  N/T  and 

(n  -  1)S2  =  E(X.  -  f.X)2/f.  =  Ef . (X./f .  -  X)2. 

i  l  i  ill  ' 

If  N  =  0,  both  numerator  and  denominator  of  V  are  zero,  and  V  is  defined  to 
have  the  value  1  for  reasons  that  will  become  clear  below. 

The  estimator  V  is  l/(n  -  1)  times  the  so-called  "index  of 
dispersion"  defined  by 

D  =  E(X.  -  f.X)2/f.X. 
l  i  i 

Note  that  D  results  from  replacing  E(Xi>  by  f .X  in  the  expression 

X2  =  E(X.  -  E(X.))2/E(X.), 

2 

and  that  D  is  also  the  X  statistic  for  testing  the  hypothesis  that  the 
observations  X.  have  Poisson  distributions.  In  this  case,  the 
conditional  distribution  of  (X  ,  X  ,  ...,  X  )  given  N  is  multinomial 


with  parameters  =  f^/T.  Hence,  under  the  Poisson  hypothesis  the  well- 

known  asymptotic  result  gives  the  conditional  distribution  of  D  given  N 
2 

as  approximately  X  with  n  -  1  degrees  of  freedom. 

This  study  was  prompted  by  an  examination  of  the  data  mentioned  in 
Sec.  I  in  which  V  was  used.  Judged  against  prevailing  beliefs  about 
values  of  p  for  actual  part  failure  distributions,  the  results  seemed 
quite  unusual — 43  of  the  hundred  values  of  V  calculated  were  greater 
than  5,  15  exceeded  15,  and  one  was  greater  than  95.  This  appeared  to 
contradict  what  was  known  about  the  properties  of  V  tor  values  of  p 
generally  considered  reasonable  and  suggested  that  these  large  values  of 
V  were  indicating  model  failures  as  well  as  larger  than  expected 
variability.  To  see  whether  model  failures  could  result  in  values  of  V 
similar  to  those  observed,  and  to  study  further  properties  of  V  under 
the  more  favorable  assumptions  described  in  Sec.  I,  V  was  examined  under 
several  sets  of  assumptions,  three  sets  of  which  will  be  discussed  here: 
that  the  number  of  failures  in  a  period  of  given  length  has  (1)  a 
Poisson  distribution  with  parameter  Xf,  (2)  a  negative  binomial 
distribution  with  parameters  r  =  Xf/(1  -  p)  and  p  =  1/p,  and  (3)  a  Poisson 
distribution  with  parameter  a  +  gf ,  a  #  0. 

Apart  from  their  relevance  to  the  part  failure  problem,  the  results 
of  (1)  are  relevant  to  inference  for  nonhomogeneous  Poisson  processes  in 
other  situations.  In  (3),  it  is  shown  that  simple  model  failure  alone 
could  not  reasonably  have  produced  the  unusual  observed  values  of  V, 
although  plots  of  failures  against  flying  hours  suggest  that  including 
an  intercept  term  in  the  negative  binomial  mean  is  desirable  in  some 
cases.  The  main  result  in  (2)  is  that  V  is  biased  low,  with  the  bias 
increasing  as  p  gets  larger  relative  to  XT.  This  implies  that  even  if 
the  negative  binomial  model  with  mean  Xf  adequately  characterizes  the 
process  generating  part  failures,  the  true  variability  is  underestimated 
on  the  average  by  V,  with  the  largest  underestimation  occurring  for  the 
largest  values  of  p. 


PROPERTIES  OF  V  UNDER  POISSON  ASSUMPTIONS 


The  properties  of  V  are  examined  here  under  Poisson  assumptions . 
Expected  values  and  variances  are  derived  for  various  values  of  n  and 
patterns  of  flying  hours,  and  Monte  Carlo  simulations  of  the 
distribution  of  V  are  displayed  and  discussed. 

The  assumptions  (henceforth  referred  to  as  "the  Poisson 
assumptions")  are  as  follows:  if  X^  is  the  number  of  demands  for  a 
given  part  in  the  ith  period,  and  f ^  the  flying  hours  for  the  ith 
period,  then  assume 

(1)  is  a  random  variable  with  a  Poisson  distribution,  having 
mean  Xf^.  That  is,  demands  are  stationary  over  periods,  with  a 
constant  expected  rate  of  demands  per  flying  hour. 

(2)  The  counts  for  different  periods  are  independent. 

Whether  the  Poisson  assumptions  hold  or  not,  the  denominator  of  V 

is  unbiased  for  X.  Under  the  Poisson  assumptions,  the  numerator  is  also 

unbiased  for  X,  and  V  has  expectation  1.  Indeed,  E(V|N)  =  1  for  all 

values  of  N  >  0.  These  results  stem  from  the  fact  that,  given  N,  the 

conditional  distribution  of  X.  is  binomial  with  parameters  N  and  p.  = 

l  A  r  x 

f./T.  It  follows  that  E(X.1N)  =  Np.  =  f.X  and 

i  i’  l 

E[(X.  -  fiX)2|N]  =  var(X.|N)  =  Np.(l  -  p^  =  f.X(l  -  p^ . 

Hence,  the  conditional  expectation  of  the  numerator  of  V  is 

E(S2!N)  =  (n  -  1)_1ZX(1  -  f^T)  =  X 

and  E(V|N)  =  E(S2/X|N)  =  1.  Taking  expectations  yields  E(S2)  =  X  and 
E(V)  =  1. 

This  generalizes  a  result  by  Rao  and  Chakravarti  (1956,  pp. 
265-266),  for  the  case  of  equal  f^.  Note  that  this  result  does  not 
require  that  the  f^  be  known  Appendix  A  contains  a  derivation  of  E(V) 
under  an  assumption  that  E(X^)  =  X^,  which  may  be  of  use  in  other 
sensitivity  analyses. 


The  above  results  do  not  hold  in  general;  i.e.,  V  is  not  usually 
unbiased  for  the  variance-to-mean  ratio.  In  particular,  it  is  biased  in 
the  negative  binomial  case,  as  will  be  seen  later. 

Both  the  conditional  and  unconditional  variances  of  V  are  of 
interest  for  assessing  the  reliability  of  V  as  an  estimator  of  p.  As  is 
shown  below,  the  conditional  variance  of  V,  given  a  nonzero  value  of  N, 
is 

var(Vi.N)  =  2(n  -  1)_1(1  -  C/N), 

where  C  =  1  +  [n2  -  TI(l/f ,)]/2(n  -  1).  If  f  =  f  =  . . .  =  f  ,  the 

i  l  z  n 

constant  C  is  equal  to  one,  and  this  is  the  largest  possible  value  of  C. 
This  follows  from  noting  that  I(l/f .)  >  I(l/f)  =  n2/T  by  Jensen's 
Inequality,  so  that 

var(VjN)  >  2(n  -  1)"1(1  -  1/N), 

with  equality  holding  if  and  only  if  the  f^  are  equal.  This  lower  bound 
should  serve  as  a  useful  approximation  in  applications  where  the  f ^  are 
approximately  equal. 

The  unconditional  variance  of  V  is  obtained  by  using  the  fact  that 

var(V)  =  var[E(V|N) J  +  E[var(V|N)], 

and  noting  that  var(V|N)  =  0  when  N  =  0  gives 

var(V)  =  2(n  -  1)_1[1  -  CE(!f  X|N  >  0)]P(N  >  0) 

=  2(n  -  l)'1!!  -  e'8  -  CH(8)] 

where  8  =  TX  and 

H(8)  =  I  "  8kc"0/k  •  k!  =  e-8  f  (eU  -  l)/u  du. 
k=l  J 

0 

Since  E(N  ^|N  >  0)  >  1/E(N|N  >  0)  =  (1  -  e  8)/8  by  Jensen's  Inequality, 


This  provides  an  upper  bound  for  var(V)  that  also  serves  as  a  good 
approximation  for  large  values  of  8.  For  a  table  of  values  of  H(8)  for 
8  between  0.01  and  20,  see  Grab  and  Savage  (1954).  The  approximation 
H(8)  -  l/(8-l)  serves  quite  well  for  8  >  5,  and  this  approximation  was 
used  in  the  work  reported  below. 

To  derive  var(V(H)  for  N  >  0,  it  is  convenient  to  set 

0  =  I(X.  -  f.X)2/f., 

l  l  i 

so  that  V  =  Q/(n  -  1)X,  and 

var(VJN)  =  var(Q|N)/(n  -  l)2^2. 

Hence,  it  suffices  to  show  that,  for  N  >  0, 

var(Q|N)  =  2(n  -  1)(1  -  C/N)X2 

where  C  is  as  above.  By  using  the  representation 

Q  =  I(X.  -  fA^/f.  -  T(X  -  X)2, 

which  implies  that 

[Q  +  T(X  -  X)2]2  =  I  U.2  +  ZU.U. 

i^j1  1  J 

2 

where  IL  =  (X^  -  f^X)  /f .,  it  is  straightforward  to  show  that 
E(Q2)  =  X[I(l/f.)  +  (1  -  2n)T]  +  X2(n2  -  1). 


E(Q2|N)  =  X[I(l/f.)  +  (1  -  2n)/T]  +  (X2  -  X/T)(n2  -  1), 

~  -2  2 
since  X  is  unbiased  for  X,  X  -  X/T  is  unbiased  for  X  ,  and  N  is  a 

2 

complete  statistic,  which  implies  that  the  unbiased  estimator  of  E(Q  ) 
that  depends  on  N  is  unique  (Lindgren,  1976,  p.  266).  To  complete  the 
derivation,  one  can  use  the  formula 

var(QlN)  =  E(Q21N)  -  [E(Q|N) ] 2 


and  the  fact  that  E(Q|N)  =  (n  -  1)X. 

To  show  how  the  values  of  var(V)  depend  on  X,  n,  and  {f_>,  four 
patterns  of  the  values  {f^}  were  chosen  so  that  the  average  was  always 
2500,  and  thus  their  sum  was  always  2500n.  These  four  patterns  are  as 
follows: 


(1)  fj  =  2500  j  =  1,  n 

(2)  fj  =  2500  j  =  2,  n  -  1 

=  1000 

f  =  4000 
n 


(3)  f ^  =  1950  +  lOOj 

f  =  1450  +  lOOj 

f.  =  1225  +  50j 

J 


for  n  =  10 
for  n  =  20 
for  n  =  50 


(4)  f  =  25  +  450j 

f  =  -  20  +  240j 
f.  =  -  50  +  lOOj 


for  n  =  10 
for  n  =  20 
for  n  =  50 


The  entries  in  the  table  below  are  var(V). 

As  this  table  shows,  the  effect  of  heterogeneity  in  flying  hours  on 
var(V)  is  more  pronounced  for  smaller  values  of  X.  Note  that  pattern 
(1)  represents  the  case  C  =  1,  and  that  the  entries  for  that  case  serve 
as  a  lower  bound  and  approximation  for  the  other  cases. 


X  =  0.01 


X  =  0.001 


f-pattem 

n  =  10 

20 

50 

n  =  10 

20 

50 

1 

0.221 

0.105 

0.041 

0.213 

0.040 

2 

0.229 

0.105 

0.041 

0.040 

3 

0.222 

0.106 

0.224 

0.115 

4 

0.242 

0.115 

0.432 

0.203 

0.116 

The  histograms  of  the  observed  values  of  V  from  selected  Monte 
Carlo  runs  appear  in  Fig.  2.  The  results  for  f^  =  2500  and  X  =  0.001 
were  essentially  the  same  as  those  for  f  =  2500  and  X  =  0.01,  so  they 
were  not  included.  One  thing  to  note  here  is  the  difference  between  the 
last  two  histograms.  Decreasing  X  from  001  to  0.001  flattens  the  mode 
and  moves  some  of  the  probability  to  each  tail- 

NEGATIVE  BINOMIAL 

The  negative  binomial  may  be  characterized  by  its  mean  and  the 
ratio  of  its  variance  to  its  mean.  This  characterization  is  of  interest 
because  the  negative  binomial  distribution  tends  to  the  Poisson  as  the 
variance-to-mean  ratio  tends  to  1,  and  thus  in  a  rough  sense  the 
variance-to-mean  ratio  describes  the  magnitude  of  the  departure  from  the 
Poisson  assumption. 

In  particular,  it  is  assumed  that  X^  has  a  negative  binomia1 
distribution  with  mean  Xf^  and  variance-to-mean  ratio  p,  so  that 

P(X^  =  k)  =  ^rj  \k  "  ^(1  -  p)Vj  k  =  0,  1,  2,  ..., 

where  r^  =  Xf^/(p  -  1),  p  =  1/p. 

The  mean  of  V  may  be  derived  by  noting  that  E(Q)  =  (n  -  l)X/p,  for 

7  7 

Q  =  (n  -  1)S“  and  S~  defined  as  in  the  introduction  to  Sec.  II.  Thus, 
E(S“)  =  X/p.  Because  X(R  +  N)/(R  +  1),  where  R  =  Ir.,  also  has 

2 

expectation  X/p,  it  follows  that  this  expression  is  equal  to  E(S~|N)  ny 
the  completeness  of  N.  Hence,  for  N  >  0, 


*  .** 


Observed  frequency  Observed  frequency  Observed  frequency 


VTMR 


Fig.  2  -  Histograms  of  V  under  Poisson  assumptions 


E (V IN)  =  E(S  /XJN)  =  (R  +  N)/(R  +  1), 


and,  as  before,  V  =  1  for  N  =  0  by  definition.  This  can  also  be  written 
as 

E(V|N)  =  1  +  (p  -  1) (N  -  1)/ (XT  +  p  -  1). 

From  either  of  these  expressions,  it  follows  immediately  that 

E(V)  =  p (R  +  P'R_1)/tR  +  1), 

which  is  less  than  p  for  p  >  1. 

This  bias  in  V  accounts  for  at  least  part  of  the  apparent 
dependence  of  variance-to-mean  ratio  on  the  mean.  As  the  mean 
increases — i.e.,  as  XT  increases — the  bias  will  decrease,  and  E(V)  will 
increase.  Thus,  even  if  all  of  the  assumptions  of  the  negative  binomial 
model  were  satisfied,  if  V  was  used  to  estimate  the  variance-to-mean 
ratio,  the  estimates  would  indicate  that  V  increases  with  the  mean. 

Although  var(V)  and  var(V|N)  are  tractable  under  the  assumptions  of 
this  subsection  using  straightforward  algebra,  the  formulas  are  very 
complicated  and  offer  little  insight. 

If  p  and  T  are  assumed  to  be  large,  the  distribution  of  V  can  be 
approximated  by  the  same  expression  when  X  is  replaced  by  X,  which  will 

•j*. 

be  called  V  .  It  can  be  shown  that  in  the  Poisson  case,  with  parameter 
Xf, 


var(V)  =  (2n  +  l  (l/f.X)]/(n  -  l)2, 

whereas,  in  the  negative  binomial  case, 

var(V*)  =  {2n/p2  +  (6q  +  p2)I(l/f .)/Xp3)]/(n  -  l)2, 

where  q  =  1  -  p.  Tnus,  when  the  f  are  large,  var(V)  in  the  negative 

2  2 

binomial  case  is  approximately  p  =  1/p  times  as  large  as  the  variance 
in  the  Poisson  case. 


Tnis  approximation  should  be  used  with  caution.  When  computed  for 
the  values  of  n,  X,  and  {f^}  used  in  the  table  of  Poisson  model 
variances  and  p  =  2,  5,  and  15,  the  approximation  was  always  greater 
than  the  exact  value,  by  factors  as  large  as  8,  with  the  error  growing 
larger  as  X  decreases,  p  increases,  or  {f^}  becomes  more  heterogeneous. 

Seme  Monte  Carlo  samples  were  drawn  from  V’s  distribution  for 
various  values  of  X  and  p,  two  patterns  of  {f.},  and  n  =  10.  Selected 
histograms  appear  in  Fig.  3(a-j).  These  figures  are  all  in  the  same 
scale  to  facilitate  making  comparisons. 

Some  generalizations  can  be  made  from  the  histograms  about  the 
effect  on  V's  distribution  of  X,  p,  and  {f^}.  First,  decreasing  X  from 
0.01  to  0.C01  shifts  some  probability  to  smaller  values  of  V,  apparently 
away  from  middle  values,  because  the  upper  tails  do  not  seem  to  differ. 
Second,  changing  from  constant  to  nonconstant  f.  also  seems  to  shift 
probability  away  from  middle  values,  but  this  is  divided  between  small 
values  and  the  upper  tail.  That  is,  changing  away  from  equal  flying 
hours  makes  more  probable  both  small  values  and  very  large  values. 

In  terms  of  the  above  trends,  the  Poisson  case  can  simply  be 
thought  of  as  the  case  with  p  =  1. 


POISSON  MODEL,  WITH  E(Xj)  =  a  +  5f. 

There  are  two  reasons  why  it  might  be  desirable  to  allow  the  mean 
of  the  failure  distribution  to  be  a  more  general  function  of  flying 
hours  than  Xf.  It  is  possible  that  some  parts  would  fail  on  planes  that 
did  not  fly  at  all,1  which  implies  that  the  function  relating  mean 
failures  to  flying  hours  can  be  approximated  by  lines  with  nonzero 
intercepts,  over  suitably  restricted  ranges  of  flying  hours.  If  flying 
hours  are  suitably  restricted,  such  approximations  could  have  either 
positive  or  negative  intercepts,  and  Fig.  4,  drawn  from  our  Air  Force 
data,  shows  that  both  of  these  possibilities  can  occur,  with  Fig.  4a 
suggest- ing  a  positive  intercept  (and  a  negative  slope),  and  Fig.  4b 
suggesting  a  negative  intercept. 


1Disuse,  especially  in  extreme  climates,  may  contribute  to  a 
significant  number  of  failures. 
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From  these  and  other  plots,  it  is  reasonable  to  postulate  a  model 
of  Poisson  failure  counts  with  mean  equal  to  a  +  Bf  for  a  not 
necessarily  zero.  But  could  observations  from  such  a  model  have 
produced  V  values  like  those  observed,  thus  providing  an  explanation  for 
their  unexpected  size?  The  results  of  this  subsection  indicate  that 
they  could  not  reasonably  be  expected  to  do  so  in  all  the  cases 
examined,  and  that  more  is  needed  to  explain  the  large  observed  V 
values. 

For  the  remainder  of  this  subsection,  the  number  of  part  failures 
in  a  given  time  period  is  assumed  to  be  a  Poisson  random  variable  with 
parameter  a  +  Bf-  Under  this  assumption,  the  formula  in  Appendix  A 
can  be  used  to  compute  E(V)  as  a  function  of  a,  B,  n,  and  {f.}. 

This  formula  is  quite  complex;  plots  of  E(V)  as  a  function  of  B, 
for  n  =  10,  for  several  values  of  a  and  for  two  sets  of  {f .} 
appear  as  Fig.  5(a-d).  Because  it  is  necessary  that  a  +  Bf.  >0, 
only  certain  values  of  B  are  permissible  when  a  <  0. 

Appendix  B  contains  a  derivation  of  var(V)  as  a  function  of  a,  B, 
n,  and  {f.}.  This  formula  is  also  quite  complex,  and  plots  of  var(V)  as 
a  function  of  B,  for  n  =  10  and  for  the  same  values  of  a  and  {f^}  as 
above,  appear  as  Fig.  6(a-d). 

A  few  trends  are  clear.  E(V)  is  larger  for  smaller  values  of  B, 
more  heterogeneous  {f^},  and  larger  |o|,  and  the  latter  effect  is  much 
more  pronounced  for  positive  a.  These  trends  also  hold  for  var(V), 
except  that  for  negative  a  and  more  heterogeneous  {f^} ,  the  size  of  B 
has  little  effect  and  that  effect  is  not  monotonic  in  B- 

Poisson  models  with  parameter  o  +  Bf^  were  fitted  to  the  100  part- 

base  combinations  in  our  data  set,  using  maximum  likelihood  estimation. 

The  results  of  these  fits  are  in  Appendix  C.  They  include  the  estimates 
*  * 

a  and  B,  their  approximate  standard  errors  and  correlation  obtained  from 
the  inverse  of  the  information  matrix,  E(V)  and  var(V)  computed  from  the 
aforementioned  formulas  using  a  and  B,  and  the  actual  observed  value  of 
V.  For  seven  part -base  combinations,  the  fitted  values  a  and  B  were 
such  that  a  +  Bf^  <  0  for  some  j;  E(V)  and  var(V)  were  not  computed  for 
those  part-base  combinations. 


fj  =  25  +  450  j 

i*  1.2.  — .  10 


fj=  1950  + 100  j 

1  =  1.2 . 10 


A  simulation  was  run  Tor  each  of  the  93  part -base  combinations 
having  all  o  +  Sf.  >0.  Observations  were  generated  for  Poisson 

^  M  A 

distributions  with  parameters  a  +  8f^.  j  =  1,2,....  n,  and  a  pseudo- 
observation  of  V  was  computed  from  them  and  compared  with  the  actual 
observed  value  of  V  300  times  for  each  admissible  part-base  combination. 
The  proportion  exceeding  the  actual  observed  value  of  V  was  recorded; 
these  proportions  appear  as  the  right -most  column  in  Appendix  C. 

The  message  of  the  simulation  and  of  E(V)  and  var(V)  is  clear. 
Almost  all  of  the  large  observed  values  of  V  are  in  the  extreme  upper 
tails  of  their  simulated  distributions,  assuming  that  failures  are 
Poisson  with  parameters  a  +  Sf^.  ^ur  ^ata  clearly  indicate  more 
variability  than  can  be  explained  by  simply  allowing  a  to  be  nonzero. 

However,  the  fitted  values  a  and  their  approximate  standard  errors 
suggest  that  a  should  not  be  assumed  to  be  zero.  Although  the  preceding 
paragraphs  cast  doubt  or.  the  appropriateness  of  a  Poisson  model,  and 
thus  cn  the  ^oprcximate  standard  errors  derived  from  that  model,  it  is 
striking  nonetheless  that  of  the  93  admissible  part-base  combinations, 

20  have  ja|  >4  s.e.(o),  28  have  |a|  >3  s.e.(a),  and  45  have  I a |  >2 
s.e.(a).  It  is  also  suggestive  that  of  the  41  admissible  part-base 
combinations  having  an  observed  V  in  excess  of  5,  17  have  |o|  >  4 
s.e.(c),  22  have  |o|),  and  29  have  jai  >  2  s.e.(a). 

This  analysis  suggests,  then,  that  some  of  the  variability  apparent 
.n  che  large  observed  values  of  V  are  a  result  of  assuming  a  =  0 
inappropriately,  but  that  the  remaining  variability  is  still  too  large 
for  a  Poisson  model  to  be  reasonable  for  all  part-base  combinations. 

CONCLUSION 

The  Poisson  model  with  parameters  Xf  is  not  adequate  for  these 
failure  data.  The  observed  values  of  V  are  too  I^rge  for  such  a 
postulated  model  to  be  credible.  Second,  these  large  values  of  V  cannot 
be  explained  reasonably  by  postulating  that  part  failures  follow  a 
Poisson  distribution  with  mean  a  +  flf  and  a  not  necessarily  zero, 
although  the  data  do  give  reason  to  believe  that  modeling  mean  failures 
this  way  is  more  appropriate  than  assuming  a  =  0. 


Taken  together,  these  two  implications  suggest  that  an  appropriate 
next  step  in  modeling  the  failure  process  would  be  to  assume  that 
failures  follow  a  negative  binomial  distribution  with  mean  a  +  flf  and 
variance-to-mean  ratio  p.  If  data  analysis  indicates  that  ct  =  0  or  p  = 
1  are  reasonable  modeling  assumptions  for  some  part-base  combination, 
then  such  a  simpler  model  may  be  appropriate  for  them. 

One  last  implication  is  that  for  cases  where  it  is  reasonable  to 
assume  that  failures  follow  a  negative  binomial  distribution  with  mean 
af  and  variance  afp,  V  may  not  be  a  desirable  estimator  for  p. 


III.  SUGGESTIONS  FOR  ESTIMATION 


NEGATIVE  BINOMIAL  WITH  MEAN  Xf 

A  simple  replacement  for  V  that  corrects  for  its  bias  has  not  been 
found.  Two  natural  alternatives  have  flaws  as  seriou  V’s.  Maximum 
likelihood  estimation,  however,  appears  to  avoid  the:.  •.  •  but  does 
not  yield  explicit  formulas  for  the  estimates.  These  ’  hree  methods  will 
now  be  examined  briefly. 

Because  E(VjN)  =  (R  +  N)/(R  -f  1)  and  the  UMVUE  of  p  based  on  N  is 
(R  +  N)/R  in  the  case  where  R  is  known,  it  may  be  desirable  to  correct  V 
by  a  multiplier  of  the  order  of  (R  +  1)/R.  One  way  to  do  this  is  to 
find  the  estimator  p  that  is  a  function  of  V  and  that,  when  V  is  set 
equal  to  its  conditional  expected  value,  reduces  to  the  desired  (R  +  N)/R 
This  estimator  is 


p  =  (N  -  1)V/(N  -  V). 

However,  it  is  easily  shown  that  the  upper  bourn!  for  V  for  given  N 

is  N(T/f  .  -  l)/(n  -  1),  where  f  .  =  min  {f.}.  This  bound  is  attained 

min  mm  j 

* 

when  X.^  =  N  for  some  j  satisfying  f =  f  .  and  X.  =  0  for  j  i  i  ; 

J~  j*  mm  j  j  j  » 

and  because  T/f  .  >  n,  the  bound  will  be  no  less  than  N.  Thus  will  be 

min  _ 

infinite  or  negative  with  positive  probability  whatever  p,  n,  X,  or  {f}  ai 
The  second  natural  suggestion  is  to  use  the  method  of  moments 
estimators  X  =  X  =  N/T  and  p  =  (IX^  -  X^F)/N,  where  F  =  Ef.^ .  However, 
it  can  be  shown  that 


Evp |N)  =  (1  -  F/T^ECVIN)  <  (1  -  n)E(V|N) 


E(p)  =  p{(l  -  p‘R'1)(l  -  F/T^R/CR  +  1)  +  P'R-1} 


<  pUR  +  p~R1)/(R  +  1)  -  R(1  -  P‘R-1)/n(R  +  1)}  <  E(V) 


for  p  >  1. 

Using  maximum  likelihood  estimates  can  provide  an  approach  to  the 


desideratum  of  the  second  paragraph  in  this  subsection.  The  log 


likelihood  function  is 


£(X,  p)  =  8  log  p  +  N  log(l  -  p)  +  £  log  |Ri  +  £i  1 J . 


where  r^  =  Xf^/(p  -  1)  and  p  =  1/p.  Reparameterize,  replacing  (X,  p)  by 


(f,  P).  For  I  =  X/(p  -  1)  and  p  =  1/p,  so  that  r^.  =  Tf^,  and  momentarily 


treat  7  as  fixed  at  I.  Then 


£(I,  p)  =  rr  log  p  +  H  log  (1  -  P)  +  K(*), 


K(I)  being  defined  in  an  obvious  way.  Differentiating  1(1,  p)  with 
respect  to  p  and  solving  for  p  yields  p  =  JT/(N  +  TT) ,  implying  that  p  = 
(N  +  TT)/2fT  =  (R  +  N)/R  for  R  =  TT.  It  remains  to  find  I,  which  can  be 
done  by  substituting  p  into  £(y,  p)  and  numerically  maximizing  the 
resulting  equation  in  I.  By  backtrans forming  and  applying  standard 
large-sample  maximum  likelihood  theory  to  the  original  log  likelihood, 
regions  may  be  derived  for  X  and  p. 


First  order  bias  corrections  and  second  order  variance 


approximations  are  available  for  these  estimators  (see,  e.g.,  Cox  and 
Hinkley,  1974,  pp.  309-310);  however,  they  involve  very  complicated 
calculations. 


NEGATIVE  BINOMIAL  WITH  MEAN  a  +  Bf 

The  log  likelihood  for  this  situation  is 


*(p,  a,  3)  =  R  log  p  +  N  log(l  -  p)  +  I  log  |Ri  "^i  1  j 


where  r^  =  (a  +  3f^)/(p  -  1)  and  p  =  1/p.  If  this  is  reparameterized  by 
setting  7  =  a/(p  -  1),  6  =  B/(p  -  1),  and  p  =  1/p,  the  reparameterized 
log  likelihood  is 


*(p,  y,  6)  =  (ny  +  6T) log  p  +  n  log  (i  -  P)  +  K(i,  6) 


K(Z,  6)  being  defined  in  an  obvious  way.  Temporarily  setting  Z  =  Z  and 
6=6  and  maximizing  in  p  yields 

p  =  (nl  +  6T)/(N  +  nZ  +  6T), 


so  that 


P  =  (N  +  nl  +  5T)/(nI  +  6T)  =  (R  +  N)/R 

for  R  =  nl  +  6T.  Thus,  maximum  likelihood  again  provides  an  approach  to 
the  desideratum  of  the  previous  subsection.  Estimates  y  and  6  can  be 
found  by  substituting  p  into  £(p,  Z,  5)  and  maximizing  numerically  in  y 
and  6.  As  before,  back  transformation  gives  estimates  for  a  and  $  and 
approximate  confidence  regions  may  be  obtained  using  standard  large 
sample  theory. 

Computing  the  method  of  moments  estimator  requires  the  solution  of 
a  system  of  three  nonlinear  equations  and  will  not  be  pursued  further. 


IV.  CONCLUSION 


The  consideration  of  competing  spare  parts  supply  systems  requires 
understanding  how  well  levels  of  parts  failures  under  wartime  conditions 
can  be  predicted.  The  ability  to  predict  levels  of  parts  failures  is 
strongly  affected  by  at  least  two  types  of  uncertainty:  about  the 
numbers  of  failures  that  will  occur  assuming  the  model  is  correct,  and 
about  the  adequacy  of  the  model  as  an  approximation  of  the  process 
generating  parts  failures. 

The  first  type  of  uncertainty  is  systematically  and  inherently 
understated  by  modeling  failures  as  Poisson  random  variables. 
Distributional  results  derived  here  make  this  quite  clear.  A  model 
that  allows  more  variability,  such  as  a  negative  binomial  model,  would 
be  more  appropriate.  The  data  also  indicate  that  prevailing  beliefs 
about  variation  within  such  negative  binomial  models  err  on  the 
optimistic  side,  and  that  inherent  variability  is  large,  for  some  parts 
at  least. 

The  second  type  of  uncertainty  can  be  accommodated  partly  by  using 
models  with  more  parameters;  this  is  the  course  taken  here  by  the 
suggestion  that  mean  failures  be  modeled  as  o  +  Df  without  assuming  a  = 
0.  Some  properties  of  possible  estimators  for  these  models  have  been 
examined  here  also.  Although  this  improvement  addresses  some  of  the 
issues  mentioned  in  Sec.  I,  it  leaves  several  issues  untouched.  To 
study  the  issues  of  stationarity  of  failure  distribution  and  the 
relationship  between  flying  hours  and  mean  failures,  more  and  better 
data  are  needed.  In  particular,  it  would  be  useful  to  have  many  years 
of  data  to  check  for  seasonality  in  failures  and  to  observe  failure 
behavior  under  a  broader  range  of  flying  hours,  especially  numbers  of 
hours  that  would  be  expected  in  wartime.  Studying  the  question  of 
independence  of  disjoint  time  intervals  awaits  the  elaboration  of 
appropriate  probabilistic  models  that  allow  for  dependence,  and  of 
suitable  estimation  methods  for  them. 


Appendix  A 


DERIVATION  CF  E(V)  UNDER  POISSON  ASSUMPTIONS 

This  appendix  derives  E(V|N)  assuming  that  X^  is  Poisson  with 
parameter  p^.  When  p^  =  Xf^,  the  Poisson  assumptions  hold  and  the 
expectations  derived  there  follow. 


rirst, 

E((Xj.  -  f ±X)2}N]  =  var(X.iS)  +  [E(X± JK)  -  fA]2 
=  Np.(l  -  p.)  +  [Np.  -  f.N/T]2 
=  Npi(l  -  pj.)  +  N^p.  -  q.)2, 

where  p^  =  p^/Ep^  and  =  f./T,  using  the  fact  that  X.jN  is  binomial 
with  parameters  N  and  p^.  It  follows  readily  that 

E(V!N)  =  (n  -  D'^NKp.  -  q.)2/qi  +  Ip.(l  -  p.)/q.} 

for  N  >  0.  Because  V  was  defined  to  be  one  for  N  =  0, 

00 

E(V)  =  P(N  =  0)  +  L  E(VjN  =  m)P(N  =  m). 
m— l 

Because  N  has  a  Poisson  distribution  with  parameter  A  =  Ip., 

E(V)  =  e"A{l  -  (n  -  l)~1I(pi  -  q.)2/q.) 

+  (n  -  i)  ^AKp.  -  qi)~/qi  +  E  p.(i  -  p^/qJ 
=  e  A  +  (n  -  1)  (A  -  e  A)I(p  -  q.)~/q.  +  Ip.  (I  -  p.)/q.}. 
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Appendix  B 

DERIVATION  OF  VAR(V/N)  AND  VAR(V)  UNDER 
POISSON  ASSUMPTIONS 


This  appendix  derives  var(V|N)  and  var(V)  under  the  assumption  that 
the  number  of  failures  in  a  given  period  has  a  Poisson  distribution  with 
parameter  y^.  The  formula  needed  for  Sec.  II  may  be  obtained  by 
substituting  y^  =  a  +  &f^. 

Using  two  facts,  that  X^jN  is  binomial  with  parameters  N  and  p^  = 
y_/Zy^,  and  that  the  joint  distribution  of  (X^,  X^  ,  N  -  X^  - 
trinomial  with  index  N  and  probabilities  p^,  ,  and  1  -  p. 

is  a  matter  of  straightforward  algebra  to  show  that 


X.  ){N  is 
J 

-  Pj  »  it 


var(ViN)  =  T^n  -  1)' 


"{4NAj  +  A2/N  +  A3), 


where 


A1  =  a5  '  a22’  A2  =  8a5  '  684  +  a3  +  4ala2  *  '  *1  ' 

2 

A3  =  -  12a5  +  6a^  -  4a3a 2  +  10a2  ,  and  a^  =  Zp^f^, 
a2  =  IPi2/fi»  a3  =  IPi/fi2>  a4  =  IPi2/fi2.  ^ 

r  3/xr  2 

a5  -  Ipi  /fi  * 

The  unconditional  variance  of  V  is  obtained  by  using 

var(V)  =  var{E(VjN)j  +  E(var(VjN)]. 

Employing  the  approximation  from  Sec.  II  and  ignoring  the  term  e  A, 
where  A  =  Zy^, 

E(var(VjN))  =  T^n  -  1)~2{4AA_  +  hj (A  -  2)  +  A.}. 

1  4  J 
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It  is  convenient  to  derive  var(E(V|N))  by  rearranging  terms  to 
avoid  difficulties  introduced  by  setting  V  =  1  when  N  =  0.  Let 


g(m)  =  (n  -  1)  1{ml(pi  -  q.)2^  +  Ip..(l  -  Pi)/q±} 


where  =  f^/T.  This  is  E(V|N  =  )  for  m  >  0.  Also,  let  G  =  E(g(N)). 

Because  N  has  a  Poisson  distribution  with  parameter  A, 

var(E(V|N))  =  £  (g(m)  -  E(V))2P(N  =  m)  +  (1  -  E(V))2e"A 

ra=l 

=  £  f g(m)  -  G]2P(N  =  m)  +  H, 

m=0 


where 


H  =  e”A[ (1  -  E(V))2  -  I  Pi(l  -  Pi)/qi  -  (G  -  E(V))2  -  (G  -  E(V))(1  -  G)] 
Thus,  var(E(VjN))  =  (n  -  1)  2A[£(pi  -  qp2/q.J2  +  aad  for  smaH  c  \ 


H  may  be  ignored- 
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Appendix  C 
Table  C.1 


ESTIMATES  AND  OTHER  QUANTITIES  FROM  FITTING  POISSON 
MODELS  WITH  MEAN  A  LINEAR  FUNCTION  OF  FLYING  HOURS 


A 

A  sv 

P(V  > 

V 

a 

B 

se(a) 

se(3) 

corr(c,g) 

E(V) 

var(V)  n  observed  V) 

(1) 

(2) 

(3) 

(4) 

(5) 

(6) 

(7) 

(8)  (9)  (10) 

.38 

-.726 

.001859 

12.700 

.002512 

-.996 

1.000 

.327 

7 

.903 

.58 

-2.851 

.006261 

7.634 

.010245 

-.998 

1.013 

.310 

7 

.777 

.68 

-2.619 

.002859 

2.472 

.000674 

-.944 

1.064 

.209 

11 

.797 

.91 

1.262 

.002112 

3.922 

.000728 

-.969 

1.015 

.174 

13 

.560 

1.00 

-4.812 

.010401 

5.772 

.007488 

-.996 

1.033 

.186 

12 

.473 

1.01 

5.595 

.001984 

4.146 

.000906 

-.953 

1.286 

.448 

9 

.597 

1.03 

4.782 

.005043 

11.173 

.001776 

-.989 

1.019 

.174 

13 

.470 

1.03 

-10.442 

.017909 

6.303 

.008252 

-.997 

1.165 

.222 

12 

.577 

1.10 

5.798 

.002566 

8.668 

.001374 

-.989 

1.044 

.184 

13 

.433 

1.22 

-1.981 

.001461 

1.554 

.000439 

-.944 

1.065 

.199 

11 

.330 

1.29 

-13.809 

.005895 

8.676 

.001431 

-.987 

1.214 

.269 

11 

.320 

1.41 

-2.980 

.002978 

2.072 

.000567 

-.928 

1.085 

.195 

12 

.220 

1.45 

-.443 

.005035 

6.859 

.008792 

-.397 

1.000 

.177 

12 

.120 

1.45 

-10.366 

.015683 

8.021 

.010391 

-.999 

■Ji. 

12 

1.50 

6.171 

-.001602 

8.091 

.010282 

-.997 

1.055 

.204 

12 

.123 

1.64 

6.708 

.001392 

7.998 

.001262 

-.991 

1.083 

.200 

13 

.107 

1.71 

25.253 

-.022562 

11.601 

.01461a. 

-.998 

1.532 

.426 

12 

.350 

1.77 

25.365 

-.000130 

10.243 

.001599 

-.991 

1.679 

.452 

13 

.373 

1.81 

1.798 

.001996 

7.533 

.001196 

-.990 

1.009 

.170 

13 

.030 

1.83 

-4.648 

.003519 

2.258 

.000622 

-.938 

1.231 

.231 

12 

.100 

1.99 

-.296 

.002393 

2.702 

.000686 

-.949 

.997 

.179 

12 

.013 

2.02 

9.758 

-.000477 

7.494 

.001275 

-.991 

1.369 

.634 

7 

.190 

2.12 

-7.822 

.003404 

6.214 

.001008 

-.987 

1.096 

.190 

13 

.027 

2.13 

-12.033 

.019514 

4.442 

.0059*0 

-.995 

1.244 

.239 

12 

.063 

2.16 

12.865 

-.001006 

3.809 

.000655 

-.981 

2.645 

1.001 

13 

.693 

2.28 

11.663 

-.000287 

3.947 

.000691 

-.975 

1.991 

.633 

13 

.320 

2.38 

-4.129 

.003868 

4.728 

.000893 

-.972 

1.051 

.177 

13 

.003 

2.41 

18.158 

-.019012 

10.257 

.913223 

- .  99S 

1 .440 

.  568 

9 

.067 

2.45 

8.537 

.003144 

12.172 

.015515 

-.997 

1.045 

.200 

12 

.003 

2.51 

15.206 

-.000205 

4.546 

.000798 

-.973 

2.186 

.708 

13 

.303 

2.54 

-7.797 

.014424 

5.963 

.007782 

-.996 

1.091 

.203 

12 

0 

2.66 

-8.396 

.004056 

3.702 

.000726 

-.967 

1.300 

.23o 

13 

.007 

2.70 

14.278 

.003035 

13.124 

.002165 

-.991 

1.155 

.30! 

10 

.013 

2.84 

1.121 

.002219 

7.523 

.001684 

-.986 

1.008 

.337 

7 

.010 

2.90 

-43.958 

.061981 

6.994 

.009586 

-.999 

* 

--- 

12 

* 

P(V  > 


V 

a 

5 

se  ( a) 

se(  B) 

corr(c,B)  E(V) 

var(V }  n  observed  V) 

(1) 

(2) 

(3) 

(A) 

(5) 

(6)  (7) 

(8)  (9)  (10) 

3 . 14 

-;.0i9 

.005797 

J  .  Ji. 

. 000874 

-.939 

1 . 146 

.218 

12 

0 

3.15 

7.797 

-.000577 

4.777 

.001487 

-.981 

1.621 

.796 

7 

.057 

3.15 

7.519 

.004659 

4.290 

.001058 

-.941 

1.326 

.327 

12 

.007 

3.15 

5.073 

.000641 

3.888 

.000701 

-.978 

1.242 

.276 

13 

.003 

3.35 

9.799 

-.000103 

3. 734 

. 000656 

-.974 

1.774 

.530 

13 

.033 

3 . 35 

63.208 

.000485 

23.886 

.004534 

-.995 

1.733 

.467 

12 

.023 

3.36 

-17.326 

.025741 

16.076 

.020651 

-1.000 

* 

X 

12 

* 

3.36 

-6.088 

.011391 

6.016 

.007813 

-.997 

1.065 

.194 

12 

.007 

3.73 

-.571 

.003320 

2.573 

.000672 

-.927 

.997 

.179 

12 

0 

3.17 

21.233 

-.001606 

4.920 

.000899 

-.973 

3.310 

1.537 

10 

.350 

3. 82 

36.181 

.000756 

25.242 

.004839 

-.996 

1.360 

.451 

9 

.007 

3.83 

-3.454 

.0G3923 

2.786 

.000765 

-.940 

1.093 

.220 

11 

0 

3.92 

-6.254 

.004965 

2.499 

.000700 

-.929 

1.311 

.258 

12 

0 

3.94 

8.479 

.000011 

4.341 

.000773 

-.962 

1.636 

.465 

13 

.010 

4.01 

-12.277 

.005804 

12.569 

.002423 

-.995 

1.090 

.210 

12 

0 

4.10 

-17.779 

.033023 

12.413 

.016036 

-.998 

1.224 

.245 

12 

0 

4.16 

-22.025 

.006090 

7.832 

.001547 

-.994 

1.550 

.343 

12 

0 

4.30 

29.394 

.004525 

8.555 

.001546 

-.971 

2.117 

.612 

13 

0 

4.42 

20.150 

.001198 

4.923 

.001203 

-.952 

3.144 

1.270 

11 

.137 

4.53 

-9.884 

.005141 

11.980 

.002308 

-.995 

1.061 

.200 

12 

0 

4.65 

-4.682 

.014296 

5.784 

.001496 

-.939 

1.043 

.193 

12 

0 

4.85 

14.496 

.000176 

10.876 

.002064 

-.995 

1.172 

.252 

12 

0 

5.07 

-124.014 

.162206 

23.615 

.030612 

-1.000 

JL 

X 

11 

JU 

5.08 

15.661 

.003134 

12.292 

.001940 

-.991 

1.183 

.237 

13 

0 

5.24 

-10.461 

.009528 

6.720 

.001281 

-.967 

1.160 

.210 

13 

0 

5.26 

.325 

.003451 

4.694 

.000877 

-.968 

1.002 

.167 

13 

0 

5.52 

-27.993 

.008101 

9.715 

.001909 

-.994 

1.616 

.368 

12 

0 

5.53 

-16.540 

.015643 

16.726 

.002678 

-.989 

1.078 

.190 

13 

0 

5.55 

14.343 

.003541 

15.628 

.002982 

-.994 

1.079 

.213 

12 

0 

5.57 

-14.772 

.012369 

14.667 

.002350 

-.989 

1.0S0 

.190 

13 

0 

5.63 

-6.142 

.002678 

6.000 

.000969 

-.989 

1.072 

.183 

13 

0 

5 . 78 

10.589 

.000478 

5.208 

.000932 

-.981 

1.628 

.445 

13 

0 

6.02 

-54.180 

.022087 

21.404 

.004143 

-.995 

1.554 

.368 

12 

0 

6.31 

-50.185 

.012110 

13.447 

.002630 

-.997 

3.186 

.708 

12 

0 

6 . 34 

-29.529 

.008996 

7.599 

.001265 

-.984 

1.762 

.350 

13 

0 

6.35 

27.142 

-  .020156 

12.328 

.015555 

-.997 

1.411 

.362 

12 

0 

6.99 

-54.733 

.012882 

11.654 

.002304 

-.997 

-U 

X 

12 

* 

7.09 

-10.435 

.005586 

8.874 

.001430 

-.988 

1.098 

.193 

13 

0 

7.28 

48.772 

-.002856 

7.248 

.001249 

-.976 

5.996 

2.487 

13 

.223 

7.38 

-104.058 

.026808 

11.370 

.001933 

-.985 

5.053 

.984 

13 

.013 

8.11 

66.438 

-.075362 

16.322 

.020428 

-.999 

* 

12 

8.25 

50.765 

-.033217 

19.743 

.025002 

-.997 

1.651 

.457 

12 

0 

8.73 

-35.961 

.052196 

13.257 

.017233 

-.999 

JL 

* 

*  n 

X  ~ 

* 

10.36 

31.387 

-.000279 

6.196 

.001084 

-.970 

3.313 

1.186 

13 

0 

0 


P(V  > 

V  a  £  se(a)  se(B)  corr(a.B)  E(V)  var(V)  n  observed  V) 


(1) 

12) 

(3) 

(4) 

(5) 

(6) 

(7) 

(8) 

(9) 

(10) 

11.11 

-48.929 

.013809 

13.114 

.002573 

-.995 

2.182 

.539 

12 

0 

11.22 

-24.905 

.015779 

6.886 

.001355 

-.957 

1.650 

.347 

13 

0 

11.34 

-27.346 

.030189 

7.181 

.001926 

-.935 

2.006 

.496 

12 

0 

13.02 

-21.252 

.013812 

3.745 

.001090 

-.930 

2.628 

.607 

12 

0 

14.09 

-55.819 

.015297 

9.574 

.001935 

-.991 

2.484 

.621 

12 

0 

14.27 

-33.828 

.019319 

21.990 

.004241 

-.994 

1.195 

.248 

12 

0 

15.05 

120.416 

-.008942 

18.086 

.002791 

-.993 

6.726 

2.782 

13 

C 

15.25 

38.252 

-.002150 

13.236 

.002063 

-.995 

2.519 

.845 

13 

0 

16.29 

43.909 

-.001820 

5.952 

.001340 

-.957 

8.441 

3.796 

12 

0 

17.94 

102.716 

-.008528 

16.106 

.002479 

-.993 

6.464 

2.729 

13 

0 

-108.728 

.027879 

12.591 

.002556 

-.993 

4.703 

1.199 

12 

0 

18.64 

-.942 

.011896 

5-819 

.001482 

-.945 

1.000 

.181 

12 

0 

19.35 

-39.378 

.010669 

9.807 

.001945 

-.994 

2.077 

.496 

12 

0 

19.84 

67.724 

-.000254 

7.823 

.001789 

-.954 

10.582 

4.690 

12 

0 

21.06 

83.637 

-.008222 

7.500 

.001233 

-.976 

13.579 

6.392 

13 

0 

22.44 

-38.009 

.024310 

29.261 

.005624 

-.996 

1.187 

.246 

12 

0 

22.46 

60.663 

-.006681 

5.570 

.000878 

-.976 

11.678 

5.604 

13 

0 

47.63 

-14.467 

.013260 

21.675 

.004163 

-.995 

1.042 

.196 

12 

0 

653.877 

-.054330 

38.079 

.005842 

-.992 

35.573 

16.230 

13 

0 

57.41 

-29.283 

.051312 

12.104 

.003102 

-.953 

1.617 

.385 

12 

0 

95.31 

-505.081 

.131377 

37.952 

.007482 

-.995 

17.436 

4.798 

12 

0 

*  indicates  that  for  that  part  and  base,  a  +  <  0 

for  some  j,  making  the  formulas  for  E(V)  and  var  (V)  inappropriate. 

Cols.  (2)  and  (3).  The  estimates  a  and  £  were  obtained 
by  maximizing  the  log  likelihood  numerically. 

Cols.  (4)  and  (5).  The  approximate  errors  were  obtained  from  the 
inverse  of  the  observed  information  matrix. 

a  ^ 

Col.  (6)  is  the  approximate  correlation  between  a  and  B. 

Cols.  (7)  and  (8).  E(V)  and  var(V)  are  computed  using  o  and  8. 

Col.  (9).  n  =  the  number  of  observations  for  that  part  and  base. 

Col.  (10)  is  the  proportion  of  300  simulated  values  of  V  greater  than 
the  observed  value  of  V,  where  the  pseudo  observations  were  generated 
from  a  Poisson  distribution  with  parameter  a  +  £f^ . 
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